	PROGRAM MOM

      use parameter_module
      use initInterface_H
      
****************************************************************************
*     INITIALIZE THE CONSTANTS
      CALL CINITAL
	CALL INITAL
****************************************************************************
*	PARAMETER INPUTS
	MUR=cmplx(1.D0,-1.d-7)
	EPSILONR=4.D0
	CONDUCTIVITY=0.D0

*******************************************************************************
*	CONSTANTS 
	X_MIN=1.D10 
	X_MAX=-1.D10
	Y_MIN=1.D10
	Y_MAX=-1.D10
*******************************************************************************
C	READ(*,*) CHOICE
	CHOICE=100

			OPEN(UNIT=150,FILE="FILL_TIME.DAT",STATUS="REPLACE")
			OPEN(UNIT=151,FILE="SOLVE_TIME.DAT",STATUS="REPLACE")

*******************************************************************************
*	INPUT STRUCTURES, ADDED BY CAO
*******************************************
	WRITE(*,*) 'ENTER THE NAME OF THE DATA TRANSFER FILE: (*.TXT)'
	READ "(A)", INFILE
C	PRINT *,'*************'
	WRITE(*,'(A)') INFILE
C	PAUSE
C	CALL READ_HEADER(EPSILONR,D,SIGMAT,SIGMAG,NUM_CELL,TANDEL,INFILE)
	CALL READ_HEADER(EPSILONR,D,SIGMAT,SIGMAG,NUM_CELL,TANDEL,
     &	ANSWER_SWEEP,FREQ_START,FREQ_END,FREQ_STEP,FREQUENCY,Z0,
     &    NUM_REGION,INFILE)

	PRINT *,'ER=',REAL(EPSILONR)
C	PRINT *,'SIGMA=',SIGMA
	PRINT *,'D=',REAL(D)
	PRINT *,'TANDEL=',REAL(TANDEL)
	PRINT *,'NUM_CELL=',NUM_CELL

	LOSSY=0
	IF(SIGMAT<1.D9) THEN
	LOSSY=1
	END IF

	GND_LOSSY=0
	IF(SIGMAG<1.D9) THEN
	GND_LOSSY=1
	END IF

	PRINT *,'LOSSY=',LOSSY
	PRINT *,'GND_LOSSY=',GND_LOSSY
	PRINT *,' '
********************************************************************************
	WRITE(*,*) 'DO YOU WANT TO DISPLAY THE CURRENT 
     &ON THE GROUND PLANE: (Y/N)?'
	READ(*,*) ANSWER_GENERATE
********************************************************************************
	IF(ANSWER_GENERATE=='Y' .OR. ANSWER_GENERATE=='Y') THEN
		WRITE(*,*)'INPUT A SINGLE FREQUENCY'
		READ(*,*) FREQUENCY
		ANS_Y_PARAM='Y'
		FREQ_STEP=0
		FREQ_MAX=FREQUENCY
	ELSE
********************************************************************************	
		WRITE(*,*) 'CURRENT SIMULATION PARAMETERS ARE:'
		IF (ANSWER_SWEEP=='Y' .OR. ANSWER_SWEEP=='Y') THEN
			WRITE(*,*) 'FREQUENY SWEEP: ON'
			WRITE(*,*) 'START FREQUENCY=',FREQ_START
			WRITE(*,*) 'END FREQUENCY=',FREQ_END
			WRITE(*,*) 'FREQUENCY STEP=',FREQ_STEP
			ANS_Y_PARAM='Y'

		FREQ_MAX=FREQ_END
		ELSE
			WRITE(*,*) 'FREQUENCY SWEEP: OFF'
			WRITE(*,*) 'SINGLE FREQUENCY=',FREQUENCY
			ANS_Y_PARAM='Y'
			FREQ_STEP=0

		FREQ_MAX=FREQUENCY
		ENDIF
		TCOUNTER_HH=1
	
		PRINT *,' '
		WRITE(*,*) 'DO YOU WANT TO CHANGE THE PARAMETERS: (Y/N)?'
		READ(*,*) ANSWER_CHANGE

********************************************************************************
		IF (ANSWER_CHANGE=='Y' .OR. ANSWER_CHANGE=='Y') THEN	
			
			WRITE(*,*) 'DO YOU WANT A FREQUENCY SWEEP: (Y/N)?'
			READ(*,*) ANSWER_SWEEP
			IF (ANSWER_SWEEP=='Y' .OR. ANSWER_SWEEP=='Y') THEN
				OPEN(UNIT=100,FILE="INPUTIMPEDANCE.DAT",STATUS="REPLACE")
				WRITE(*,*)'INPUT THE START FREQUENCY (HZ)'
				READ(*,*) FREQ_START
				WRITE(*,*)'INPUT THE END FREQUENCY (HZ)'
				READ(*,*) FREQ_END
				WRITE(*,*)'INPUT THE NUMBER OF STEPS'
				READ(*,*) FREQ_STEP
				ANS_Y_PARAM='Y'
				FREQUENCY=FREQ_START
			FREQ_MAX=FREQ_END
			ELSE
C				FREQUENCY=8.D9
				WRITE(*,*)'INPUT A SINGLE FREQUENCY'
				READ(*,*) FREQUENCY
				WRITE(*,*)'DO YOU WANT S-PARAMETERS? (Y/N)'
				READ(*,*) ANS_Y_PARAM
				FREQ_STEP=0

			FREQ_MAX=FREQUENCY
			ENDIF
	
		END IF

	END IF

	TCOUNTER_HH=1
		
		IF (ANS_Y_PARAM=='Y' .OR. ANS_Y_PARAM=='Y') THEN
C			OPEN(UNIT=111,FILE="OUTPUT_Y11.DAT",STATUS="REPLACE")
C			OPEN(UNIT=112,FILE="OUTPUT_Y12.DAT",STATUS="REPLACE")
C			OPEN(UNIT=113,FILE="OUTPUT_Y21.DAT",STATUS="REPLACE")
C			OPEN(UNIT=114,FILE="OUTPUT_Y22.DAT",STATUS="REPLACE")
			OPEN(UNIT=115,FILE="OUTPUT_S11.DAT",STATUS="REPLACE")
			OPEN(UNIT=116,FILE="OUTPUT_S21.DAT",STATUS="REPLACE")
			OPEN(UNIT=117,FILE="OUTPUT_S12.DAT",STATUS="REPLACE")
			OPEN(UNIT=118,FILE="OUTPUT_S22.DAT",STATUS="REPLACE")
C			OPEN(UNIT=119,FILE="OUTPUT_S13.DAT",STATUS="REPLACE")
C			OPEN(UNIT=120,FILE="OUTPUT_S23.DAT",STATUS="REPLACE")
C			OPEN(UNIT=121,FILE="OUTPUT_S14.DAT",STATUS="REPLACE")
C			OPEN(UNIT=122,FILE="OUTPUT_S24.DAT",STATUS="REPLACE")
C			OPEN(UNIT=123,FILE="OUTPUT_S31.DAT",STATUS="REPLACE")
C			OPEN(UNIT=124,FILE="OUTPUT_S41.DAT",STATUS="REPLACE")
C			OPEN(UNIT=125,FILE="OUTPUT_S32.DAT",STATUS="REPLACE")
C			OPEN(UNIT=126,FILE="OUTPUT_S42.DAT",STATUS="REPLACE")
C			OPEN(UNIT=127,FILE="OUTPUT_S33.DAT",STATUS="REPLACE")
C			OPEN(UNIT=128,FILE="OUTPUT_S43.DAT",STATUS="REPLACE")
C			OPEN(UNIT=129,FILE="OUTPUT_S34.DAT",STATUS="REPLACE")
C			OPEN(UNIT=130,FILE="OUTPUT_S44.DAT",STATUS="REPLACE")
			WRITE(*,*) 'HOW MANY PORTS DOES THE SYSTEM HAVE?'
			READ(*,*) NUM_PORT
*       TOUCHSTONE OUTPUT
			IF(NUM_PORT==1) THEN
			OPEN(UNIT=155,FILE="OUTPUT_S.S1P",STATUS="REPLACE")
			ELSE IF(NUM_PORT==2) THEN
			OPEN(UNIT=155,FILE="OUTPUT_S.S2P",STATUS="REPLACE")
			ENDIF
			WRITE(155,*) '# HZ S RI R',Z0
		ENDIF

********************************************************************************
	LOSSY=0
	IF(SIGMAT<1.D9) THEN    ! SIGMAT: SIGMA TRACE
	LOSSY=1
	END IF

	GND_LOSSY=0
	IF(SIGMAG<1.D9) THEN
	GND_LOSSY=1
	END IF
C	PRINT *,'LOSSY=',LOSSY
C	PRINT *,'GND_LOSSY=',GND_LOSSY
	ALLOCATE(GLOBAL_HEIGHT(0:NUM_REGION))
	ALLOCATE(GLOBAL_PERMITTIVITY(0:NUM_REGION))
	ALLOCATE(GLOBAL_TANDEL(0:NUM_REGION))
	ALLOCATE(REGION(1:NUM_CELL,1:4))         !??
	ALLOCATE(GLOBAL_K(0:NUM_REGION))


	CALL READ_STRUCT(NUM_CELL,CELL_Z,PORT1,PORT2,
     &	PORT11,PORT22,PC1,PC2,PC11,PC22,LOAD,SOURSE,ZC,SC,
     &	LDS,SS,Z_VALUE,S_VALUE,REGION,NUM_REGION,
     &    GLOBAL_HEIGHT,GLOBAL_PERMITTIVITY,GLOBAL_TANDEL,INFILE)	

	SEGS_W=PC1

C	DO I=1,NUM_CELL
C	CELL_Z(I)%SECTION=1
C	CELL_Z(I)%LAYER=1
C	END DO
C   To check if alpha is a special case
	DO I=1,NUM_CELL
	  DO J=1,4
	    IF(ABS(CELL_Z(I)%ALPHA(J)-PI/2.D0)<1.D-6) THEN
	      CELL_Z(I)%ALPHA(J)=PI/2.D0
	    ELSE IF(ABS(CELL_Z(I)%ALPHA(J)+PI/2.D0)<1.D-6) THEN
	      CELL_Z(I)%ALPHA(J)=-PI/2.D0
	    ELSE IF(ABS(CELL_Z(I)%ALPHA(J)-PI)<1.D-6) THEN
	      CELL_Z(I)%ALPHA(J)=PI
	    ELSE IF(ABS(CELL_Z(I)%ALPHA(J)+PI)<1.D-6) THEN
	      CELL_Z(I)%ALPHA(J)=-PI
	    ELSE IF(ABS(CELL_Z(I)%ALPHA(J)-PI*3.D0/2.D0)<1.D-6) THEN
	      CELL_Z(I)%ALPHA(J)=PI*3.D0/2.D0
	    ELSE IF(ABS(CELL_Z(I)%ALPHA(J)+PI*3.D0/2.D0)<1.D-6) THEN
	      CELL_Z(I)%ALPHA(J)=-PI*3.D0/2.D0
	    END IF
	  END DO
	END DO

	OPEN(2,FILE='INPUTCONFIRM.TXT',STATUS='UNKNOWN')

	WRITE(2,'("ER=",F10.6,X,"SIG=",F18.6,X,"D=",F12.6,X,"NUM_CELL="
     &	,I4)'),EPSILONR,SIGMAT,D,NUM_CELL
     	DO I=1,NUM_CELL
	WRITE(2,'(I3,I3,I3,I3,4I3,4F21.6,4F21.6,4F21.6,4F21.6,F21.6,4F21.6,
     &	A4,A4,I2,I2)'),CELL_Z(I)%SECTION,CELL_Z(I)%LAYER,
     &    CELL_Z(I)%FLAG,
     &	CELL_Z(I)%NUM_SUBCELL,CELL_Z(I)%CURR_FLAG,
     &	CELL_Z(I)%L_OR_TANPHI,CELL_Z(I)%B,CELL_Z(I)%X,CELL_Z(I)%Y,
     &	CELL_Z(I)%Z,CELL_Z(I)%ALPHA,CELL_Z(I)%CODE,CELL_Z(I)%LR_CODE,
     &	LEN(CELL_Z(I)%CODE),LEN(CELL_Z%LR_CODE)
C	PRINT '(I3,/)',I
C	PRINT '(I3,I3,F20.6,/)',CELL_Z(I)%FLAG,CELL_Z(I)%NUM_SUBCELL,
C     &	CELL_Z(I)%B(1)	
	ENDDO

	WRITE(2,'("ER=",F10.3,X,"D=",F10.6)'),EPSILONR,D
	DO I=1,PC1
	WRITE(2,'("PORT_1",I2,"=",I6)'),I,PORT1(I)
	ENDDO
	DO I=1,PC11
	WRITE(2,'("PORT_11",I2,"=",I6)'),I,PORT11(I)
	ENDDO
	DO I=1,PC2
	WRITE(2,'("PORT_2",I2,"=",I6)'),I,PORT2(I)
	ENDDO
	DO I=1,PC22
	WRITE(2,'("PORT_22",I2,"=",I6)'),I,PORT22(I)
	ENDDO

	DO I=1,LDS
	DO J=1,ZC(I)
	WRITE(2,'("LOAD",I2,I2,"=",I6," VALUE:")'),I,J,LOAD(I,J)
     	WRITE(2,*),Z_VALUE(I)
	END DO
	END DO

	DO I=1,SS
	DO J=1,SC(I)
	WRITE(2,'("SOURSE",I2,I2,"=",I6)'),I,J,SOURSE(I,J)
	END DO
	END DO


	CLOSE(2)

	

***    PRECOMPUTE RHO_MAX TO INDICATE WHETHER ZZ_1 NEED TO BE COMPUTED-------------
      ALLOCATE ( INDICATION (1:NUM_CELL) )

	MU=MUO*MUR

C	TANDEL=CONDUCTIVITY/W/EPSILONO/EPSILONR

	TANDEL=0.0001D0

	EPSILON=EPSILONO*EPSILONR*(1.D0-IM*TANDEL)

*********     CALCULATE THE MAXIMUM RADIUS WHERE WE NEED TO CALCULATE HIGHER ORDER MODES.
*********     ALL HORIZONTAL IMPEDANCE ELEMENTS CAN BE SET TO ZERO IF RHO>RHO_MAX
	RHO_MAX=ABS(LOG(RELERR_RHOMN_HH)/DBLE(SQRT((PI/D)**2
     &            	-(TWOPI*FREQ_MAX)**2* MU*EPSILON) ))

      DO I=1, NUM_CELL
      INDI_HEAD => INDICATION (I)
	NULLIFY (INDI_HEAD%NEXT)
	END DO



	DO I=1, NUM_CELL
C	DO I=4,4
*		GET TEST FUNCTION DATA
C	  PAUSE
	  Z_M=CELL_Z(I)%Z
	  SEC_M=CELL_Z(I)%SECTION
	  LAYER_M=CELL_Z(I)%LAYER
	  FLAG_M=CELL_Z(I)%FLAG
	  SUBCELL_M=CELL_Z(I)%NUM_SUBCELL
C        PP=1
	  DO J= I, NUM_CELL
C	  DO J=5,5
*		GET EXPANSION FUNCTION DATA	
		Z_N=CELL_Z(J)%Z
		SEC_N=CELL_Z(J)%SECTION
		LAYER_N=CELL_Z(J)%LAYER
		FLAG_N=CELL_Z(J)%FLAG
		SUBCELL_N=CELL_Z(J)%NUM_SUBCELL

	    CODESIZE=0

		DO III=1,16
		  CODE_MN(III)=' '
		  REAC_TYPE(III)=' '

		  TEST_LTANPHI(III)=0.D0
		  EXP_LTANPHI(III)=0.D0
		  TEST_CURR(III)=0
		  EXP_CURR(III)=0
		  TEST_X(III)=0.D0
		  EXP_X(III)=0.D0
		  TEST_Y(III)=0.D0
		  EXP_Y(III)=0.D0
		  TEST_B(III)=0.D0
		  EXP_B(III)=0.D0
		  TEST_ALPHA(III)=0.D0
		  EXP_ALPHA(III)=0.D0	
		ENDDO			


		CALL DECODE_CELL(CELL_Z(I),CELL_Z(J),CODESIZE,CODE_MN,
     &					 TEST_LTANPHI,EXP_LTANPHI,TEST_X,EXP_X,
     &					 TEST_Y,EXP_Y,TEST_B,EXP_B,TEST_ALPHA,
     &				     EXP_ALPHA,TEST_CURR,EXP_CURR,REAC_TYPE,
     &					 SUBREAC_TYPE)

************************************************************************************ 
****  RIGHT NOW, WE ONLY CONSIDER THE FIRST OF TEST_LTANPHI
************************************************************************************ 
		DO III=1,CODESIZE	  
		  X_M=TEST_X(III)

		  IF (X_M < X_MIN) THEN
			X_MIN = X_M
		  END IF

		  IF (X_M > X_MAX) THEN
			X_MAX = X_M
		  END IF

		  Y_M=TEST_Y(III)

		  IF (Y_M < Y_MIN) THEN
			Y_MIN = Y_M
		  END IF
		
		  IF (Y_M > Y_MAX) THEN
			Y_MAX = Y_M
		  END IF			

		  B_M=TEST_B(III)
		  L_M=TEST_LTANPHI(III)
		  ALPHA_M=TEST_ALPHA(III)

		  X_N=EXP_X(III)
		  Y_N=EXP_Y(III)
			
		  B_N=EXP_B(III)
		  L_N=EXP_LTANPHI(III)
		  ALPHA_N=EXP_ALPHA(III)
		  
		  CURR_FLAG_M=TEST_CURR(III)
		  CURR_FLAG_N=EXP_CURR(III)
*  GET CURRENT FLAG MULTIPLYER
		  CURR_FLAG_MN=CURR_FLAG_M*CURR_FLAG_N

		  ALPHA_MN=ALPHA_M-ALPHA_N
		  COS_ALPHA_MN=COS(ALPHA_MN)
		  SIN_ALPHA_MN=SIN(ALPHA_MN)

		  IF (ALPHA_MN .EQ. 0.D0) THEN
			ALPHA_MN_EQ0=.TRUE.
		  ELSE
			ALPHA_MN_EQ0=.FALSE.
	      END IF

** COMPUTE THE AVERAGE DISTANCE BETWEEN THE TWO REACTION CELLS
		  X_MN=X_M*COS_ALPHA_MN-Y_M*SIN_ALPHA_MN-X_N
		  Y_MN=X_M*SIN_ALPHA_MN+Y_M*COS_ALPHA_MN-Y_N
		  
		  RHO_CENTER=SQRT(X_MN*X_MN+Y_MN*Y_MN)
		  
		  RHO_MIN=RHO_CENTER-B_M*7.D0
		  
		  THETA_MN_CENTER=ATAN2(-Y_MN,X_MN)

************************************************************************************ 
***	TEST TO SEE IF THE COUPLING IS NEGLIGIBLE
************ IN THE FOLLOWING SITUATIONS, ZZ_1 NEEDS TO BE CALCULATED
	      IF (RHO_CENTER <= 2.1D0 * MAX(L_N, L_M) ) THEN
			INDI_HEAD => INDICATION(I)
			CALL INDICATION_FILL(INDICATION(I), J, III)

			 ELSEIF (RHO_MIN .LE. RHO_MAX*0.8D0) THEN
			 INDI_HEAD => INDICATION(I)
			 CALL INDICATION_FILL(INDICATION(I), J, III)
	      END IF
*****************************************************************
          END DO  !!!!!!!!! END III
	 END DO     !!!!!!!!! END J
	END DO		!!!!!!!!! END I

****-------------------------------------------------------
	PRINT '("NOW COMPUTING AND SOLVING THE REACTION MATRIX",/)'
C	PAUSE

*** START THE FREQUENCY SWEEP LOOP*************
      TIME2 = CPSEC()

	DO KK=1,FREQ_STEP+1

      ALLOCATE (AA_CG(1:NUM_CELL))
	ALLOCATE (R_CG(1:NUM_CELL))
	ALLOCATE (TMP_CG(1:NUM_CELL))
	ALLOCATE (P_CG(1:NUM_CELL))
	ALLOCATE (Q_CG(1:NUM_CELL))
	ALLOCATE  (Z_CG(1:NUM_CELL))

	LAMBDA=3.D8/SQRT(EPSILONR)/FREQUENCY

** READ IN THE SIMULATION STRUCTURE
** THE STRUCTURE MAY VARY WITH FREQUENCY

	IF (CHOICE==1) THEN
*	HORIZONTAL DIPOLE TEST STRUCTURE
		CALL STRUCT_HORIZONTAL_DIPOLE(FREQUENCY,EPSILONR,LAMBDA,CELL_Z,
     &				D,NUM_CELL,MAX_CELLS)
	ELSEIF (CHOICE==2) THEN
*		COUPLED AND CROSSOVER MULIPLE LINES
		CALL STRUCT_MULTI_LINE(FREQUENCY,EPSILONR,LAMBDA,CELL_Z,D,
     &		NUM_CELL,MAX_CELLS)
	ELSEIF (CHOICE==3) THEN
*		TEST STRUCTURE FOR IMPROVING EFFICIENCY
		CALL STRUCT_TEST(FREQUENCY,EPSILONR,LAMBDA,CELL_Z,D,
     &		NUM_CELL,MAX_CELLS)
	ELSEIF (CHOICE==4) THEN
*		TEST STRUCTURE WITH HALF CELL
		CALL STRUCT_TEST_HALFCELL(FREQUENCY,EPSILONR,LAMBDA,CELL_Z,D,
     &		NUM_CELL,MAX_CELLS)
	ELSEIF (CHOICE==5) THEN
*	VERTICAL DIPOLE TEST STRUCTURE
		CALL STRUCT_VER_DIPOLE_WANG(FREQUENCY,EPSILONR,LAMBDA,
     &	                    CELL_Z,D, NUM_CELL,MAX_CELLS)
	ELSEIF (CHOICE==6) THEN
*	SINGLE TRANSMISSION LINE WITH VIAS 
		CALL  STRUCT_TRANS_LINE_WANG(FREQUENCY,EPSILONR,LAMBDA,CELL_Z
     &        	,D, NUM_CELL,MAX_CELLS)
	ELSEIF (CHOICE==7) THEN
*		DD TEST STRUCTURE
		CALL STRUCT_DD(FREQUENCY,EPSILONR,LAMBDA,CELL_Z,
     &				D,NUM_CELL,MAX_CELLS)

	ENDIF
*******************************************************************************
*	CALCULATIONS BASED ON INPUT DATA

	W=TWOPI*FREQUENCY
	MU=MUO*MUR

	TANDEL=0.0001D0

	EPSILON=EPSILONO*EPSILONR*(1.D0-IM*TANDEL)

	K=W*SQRT(MU*EPSILON)

*     CALCULATE THE MAXIMUM RADIUS WHERE WE NEED TO CALCULATE HIGHER ORDER MODES.
*     ALL HORIZONTAL IMPEDANCE ELEMENTS CAN BE SET TO ZERO IF RHO>RHO_MAX
	RHO_MAX=ABS(LOG(RELERR_RHOMN_HH)/DBLE(SQRT((PI/D)**2-K*K)))

******************************************************************************
*	PRINT STRUCTURE DATA TO FILE  (KK for the frequency step)
	IF (KK==1) THEN
		OPEN(UNIT=110,FILE="STRUCTURE.TXT",STATUS="REPLACE")
		WRITE(110,*) 'LAMBDA=',LAMBDA	
		WRITE(110,*) 'EPSILONR=',EPSILONR
		WRITE(110,*) 'MUR=',MUR
		WRITE(110,*) 'LOSS TANGENT=',TANDEL
		WRITE(110,*) 'D=',D
		WRITE(110,*) 'NUMBER OF CELLS: ',NUM_CELL
		WRITE(110,*) ' '

		DO I=1,NUM_CELL
			WRITE(110,*)'CELL #',I
			WRITE(110,*)'Z=',CELL_Z(I)%Z
			WRITE(110,*)'SECTION=',CELL_Z(I)%SECTION
			WRITE(110,*)'LAYER=',CELL_Z(I)%LAYER
			WRITE(110,*)'FLAG=',CELL_Z(I)%FLAG
			WRITE(110,*)'CODE=',CELL_Z(I)%CODE
			WRITE(110,*)'LEFT-RIGHT CODE=',CELL_Z(I)%LR_CODE
			WRITE(110,*)'CURRENT FLAG(1ST)=',CELL_Z(I)%CURR_FLAG(1)
			WRITE(110,*)'CURRENT FLAG(2ND)=',CELL_Z(I)%CURR_FLAG(2)
			WRITE(110,*)'CURRENT FLAG(3RD)=',CELL_Z(I)%CURR_FLAG(3)
			WRITE(110,*)'CURRENT FLAG(4TH)=',CELL_Z(I)%CURR_FLAG(4)
			WRITE(110,*)'# OF SUBCELL:',CELL_Z(I)%NUM_SUBCELL
			DO J=1,CELL_Z(I)%NUM_SUBCELL
				WRITE(110,*)'SUBCELL #:',J
				WRITE(110,*)'X=',CELL_Z(I)%X(J)
				WRITE(110,*)'Y=',CELL_Z(I)%Y(J)
				WRITE(110,*)'L=',CELL_Z(I)%L_OR_TANPHI(J)
				WRITE(110,*)'B=',CELL_Z(I)%B(J)
				WRITE(110,*)'ALPHA=',CELL_Z(I)%ALPHA(J)*180/PI
				WRITE(110,*)' '
			END DO 
		END DO 
		
		CLOSE(110)		
	ENDIF
*	STOP '***STRUCTURE FILE DOWNLOADED***'
C	WRITE(*,*) 'FREQUENCY=',FREQUENCY
*******************************************************************************
*	RUN MOM CODE

      DO KM=1,MAX_PLATES
	DO LAYER_M=1,MAX_LAYER
	 DO LAYER_N=1,MAX_LAYER
*     INITIALIZE THE GREENS FUNCTION ARRAY COUNTERS FOR THE VARIOUS LAYERS
        NUM_GREEN_HH(LAYER_M,LAYER_N, KM)=0
	 ENDDO
	ENDDO
	ENDDO

***!!!!!!!!!!!!!!!!! INITIALIZE ALL ELEMENTS OF SEVERAL LOOK_UP_TABLE
	DO KM=1, MAX_PLATES
	DO SEC_M=1,100
	 DO SEC_N=1,100
      
	HEAD => Z_LOOKUP_RRHH(SEC_M, SEC_N, KM)
	NULLIFY (HEAD%NEXT)	

	HEAD => Z_LOOKUP_RFHH(SEC_M, SEC_N, KM)
	NULLIFY (HEAD%NEXT)	
		   

	HEAD => Z_LOOKUP_HH_0(SEC_M, SEC_N, KM)
	NULLIFY (HEAD%NEXT)

	HEAD => Z_LOOKUP_HH_PI(SEC_M, SEC_N, KM)
	NULLIFY (HEAD%NEXT)

	HEAD => Z_LOOKUP_HH_90_PLUS(SEC_M, SEC_N, KM)
	NULLIFY (HEAD%NEXT)	

	HEAD => Z_LOOKUP_HH_90_MINUS(SEC_M, SEC_N, KM)
	NULLIFY (HEAD%NEXT)	
*-----------
	HEAD2 => Z_LOOKUP_KK_0(SEC_M, SEC_N, KM)
	NULLIFY (HEAD2%NEXT)

	HEAD2 => Z_LOOKUP_KK_PI(SEC_M, SEC_N, KM)
	NULLIFY (HEAD2%NEXT)

	HEAD2 => Z_LOOKUP_KK_90_PLUS(SEC_M, SEC_N, KM)
	NULLIFY (HEAD2%NEXT)	

	HEAD2 => Z_LOOKUP_KK_90_MINUS(SEC_M, SEC_N, KM)
	NULLIFY (HEAD2%NEXT)
	
	HEAD1 => Z_LOOKUP_KH_0(SEC_M, SEC_N, KM)
	NULLIFY (HEAD1%NEXT)

	HEAD1 => Z_LOOKUP_KH_PI(SEC_M, SEC_N, KM)
	NULLIFY (HEAD1%NEXT)

	HEAD1 => Z_LOOKUP_KH_90_PLUS(SEC_M, SEC_N, KM)
	NULLIFY (HEAD1%NEXT)	

	HEAD1 => Z_LOOKUP_KH_90_MINUS(SEC_M, SEC_N, KM)
	NULLIFY (HEAD1%NEXT)	
*-----------
	 ENDDO
	ENDDO
	ENDDO


      DO I=1, NUM_CELL
      AA_CG_HEAD => AA_CG(I)
	NULLIFY(AA_CG_HEAD%NEXT)
      END DO

** CURRENTLY ONLY TAKE CARE THE HOMOGENEOUS GREEN'S FUNCTION.

* FIND KRHO FOR P=0 MODE
	NUM_KRHO=0

	DO KM=1, NUM_REGION
	   TANDEL          = GLOBAL_TANDEL(KM)
	   EPSILONR        = GLOBAL_PERMITTIVITY(KM)
         EPSILON         = EPSILONO*EPSILONR*(1.D0-IM*TANDEL)
         GLOBAL_K(KM)    = W*SQRT(MU*EPSILON)
         STORE_KRHO(0,KM)  = GLOBAL_K(KM)
	END DO

	DO I=1, NUM_CELL

		Z_M=CELL_Z(I)%Z - GLOBAL_HEIGHT(REGION(I,1)-1)
		LAYER_M=CELL_Z(I)%LAYER

		DO J= 1, NUM_CELL 

			Z_N=CELL_Z(J)%Z - GLOBAL_HEIGHT(REGION(J,1)-1)
			LAYER_N=CELL_Z(J)%LAYER

			IF (REGION(I,1)==REGION(J,1)) THEN

	  REGION_M = REGION(I,1)
	    TANDEL = GLOBAL_TANDEL(REGION_M)
	  EPSILONR = GLOBAL_PERMITTIVITY(REGION_M)
         EPSILON = EPSILONO*EPSILONR*(1.D0-IM*TANDEL)
               K = GLOBAL_K(REGION_M)
	         D = GLOBAL_HEIGHT(REGION_M)- GLOBAL_HEIGHT(REGION_M-1)
		
			CALL COMPUTE_GREEN(Z_M,Z_N,K,D,
     &				NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &				NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &				FREQUENCY,SIGMAG,GND_LOSSY,EPSILON,MU,
     &				STORE_KRHO_TM(0,REGION_M),
     &                STORE_KRHO_TE(0,REGION_M))

			END IF

		END DO
	END DO
**************************************************************************
*	START TIMER FOR MATRIX FILL TIME
      TIME0 = CPSEC()

**** I IS THE LOOP FOR TESTING CELL
**** J IS THE LOOP FOR EXPANSION CELL

C      GO TO 1000	
      OPEN(UNIT=60,FILE="REACTION.DAT", STATUS="REPLACE")

      INDEX=1
      
	DO I=1, NUM_CELL
C      DO I=1,40
*		GET TEST FUNCTION DATA
C	  PAUSE
	  Z_M=CELL_Z(I)%Z
	  SEC_M=CELL_Z(I)%SECTION
	  LAYER_M=CELL_Z(I)%LAYER
	  FLAG_M=CELL_Z(I)%FLAG
	  SUBCELL_M=CELL_Z(I)%NUM_SUBCELL	
	  
*	  WE ONLY NEED TO FILL THE UPPER HALF OF THE Z-MATRIX BECAUSE OF RECIPROCITY

      INDI_HEAD  => INDICATION(I)
	INDI_TRAIL => INDI_HEAD
	INDI_TMP   => INDI_HEAD%NEXT
	  
	J= INDI_TMP%J
	ZZ_1=0

	  DO    J=I, Num_cell 
*---------=======+++++++++++++++++++===========------------       
         ZZ_1=0
	
		Z_N=CELL_Z(J)%Z
		SEC_N=CELL_Z(J)%SECTION
		LAYER_N=CELL_Z(J)%LAYER
		FLAG_N=CELL_Z(J)%FLAG
		SUBCELL_N=CELL_Z(J)%NUM_SUBCELL
*****	CLEAR DECODE_CELL OUTPUTS
		CODESIZE=0
		DO III=1,16
		  CODE_MN(III)=' '
		  REAC_TYPE(III)=' '
		  TEST_LTANPHI(III)=0.D0
		  EXP_LTANPHI(III)=0.D0
		  TEST_CURR(III)=0
		  EXP_CURR(III)=0
		  TEST_X(III)=0.D0
		  EXP_X(III)=0.D0
		  TEST_Y(III)=0.D0
		  EXP_Y(III)=0.D0
		  TEST_B(III)=0.D0
		  EXP_B(III)=0.D0
		  TEST_ALPHA(III)=0.D0
		  EXP_ALPHA(III)=0.D0	
		ENDDO			

		CALL DECODE_CELL(CELL_Z(I),CELL_Z(J),CODESIZE,CODE_MN,
     &					 TEST_LTANPHI,EXP_LTANPHI,TEST_X,EXP_X,
     &					 TEST_Y,EXP_Y,TEST_B,EXP_B,TEST_ALPHA,
     &				     EXP_ALPHA,TEST_CURR,EXP_CURR,REAC_TYPE,
     &					 SUBREAC_TYPE)

*** CLEAR REACTION ELEMENT
		  III = INDI_TMP %III
	    INDI_TRAIL => INDI_TMP
999	    INDI_TMP   => INDI_TMP%NEXT

*** NEED TO ADD IN NEW TERMS TO THE NEW VERSIONS LOOP

        DO KM=1, SUBCELL_M
	  DO KN=1, SUBCELL_N

	     III = (KM-1)*SUBCELL_N +KN
	     REGION_M = REGION(I, KM)
	     REGION_N = REGION(J, KN)
*          IF TWO CELLS ARE IN THE SAME REGION(PLATE), THEN PROCESSING.
*          IF two cells are not in the same region, not calculated.
	     IF (REGION_M == REGION_N) THEN

             Z_M = CELL_Z(I)%Z - GLOBAL_HEIGHT(REGION_M-1)
             Z_N = CELL_Z(J)%Z - GLOBAL_HEIGHT(REGION_M-1)
	 
	    TANDEL = GLOBAL_TANDEL(REGION_M)
	  EPSILONR = GLOBAL_PERMITTIVITY(REGION_M)
         EPSILON = EPSILONO*EPSILONR*(1.D0-IM*TANDEL)
               K = GLOBAL_K(REGION_M)
	         D = GLOBAL_HEIGHT(REGION_M)- GLOBAL_HEIGHT(REGION_M-1)

         RHO_MAX = ABS(LOG(RELERR_RHOMN_HH)/DBLE(SQRT((PI/D)**2-K*K)))
*** SET TRANSFER VALUES===============
		  X_M=TEST_X(III)
		  Y_M=TEST_Y(III)
		  B_M=TEST_B(III)
		  L_M=TEST_LTANPHI(III)
		  ALPHA_M=TEST_ALPHA(III)

		  X_N=EXP_X(III)
		  Y_N=EXP_Y(III)
		  B_N=EXP_B(III)
		  L_N=EXP_LTANPHI(III)
		  ALPHA_N=EXP_ALPHA(III)
		  
		  CURR_FLAG_M=TEST_CURR(III)
		  CURR_FLAG_N=EXP_CURR(III)
*  GET CURRENT FLAG MULTIPLYER
		  CURR_FLAG_MN=CURR_FLAG_M*CURR_FLAG_N

		  ALPHA_MN=ALPHA_M-ALPHA_N
		  COS_ALPHA_MN=COS(ALPHA_MN)
		  SIN_ALPHA_MN=SIN(ALPHA_MN)

		  IF (ALPHA_MN .EQ. 0.D0) THEN
			ALPHA_MN_EQ0=.TRUE.
		  ELSE
			ALPHA_MN_EQ0=.FALSE.
	    END IF

** COMPUTE THE AVERAGE DISTANCE BETWEEN THE TWO REACTION CELLS
		  X_MN=X_M*COS_ALPHA_MN-Y_M*SIN_ALPHA_MN-X_N
		  Y_MN=X_M*SIN_ALPHA_MN+Y_M*COS_ALPHA_MN-Y_N
		  
		  RHO_CENTER=SQRT(X_MN*X_MN+Y_MN*Y_MN)
		  
		  IF ((X_MN.EQ.0.D0).AND.(Y_MN.EQ.0.D0)) THEN
			THETA_MN_CENTER=-PI/2.D0
		  ELSE
			THETA_MN_CENTER=ATAN2(-Y_MN,X_MN)
		  END IF
		  
		  RHO_MIN=RHO_CENTER-B_M*7.D0

***	TEST TO SEE IF THE COUPLING IS NEGLIGIBLE
CXX		  IF (RHO_MIN .LE. RHO_MAX*4.D0) THEN
		  		  
		  SELECT CASE (REAC_TYPE(III)(2:2))
		    CASE ('H')
*			  CASES _H__
			  SELECT CASE (REAC_TYPE(III)(4:4))
			    CASE ('H')
*			      CASES _H_H
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('B')
*					  CASES BH_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE BHFH
		   				  ZZ_1=ZZ_1+NI_BHFH
     &	(TEST_X,TEST_Y,Z_M,TEST_LTANPHI,TEST_B,TEST_ALPHA,EXP_X,EXP_Y,
     &	Z_N,EXP_LTANPHI,EXP_B,EXP_ALPHA,SUBREAC_TYPE,K,D,W,EPSILON,
     &	RHO_MAX)			  	
						CASE ('B')
*						  CASE BHBH
		   				  ZZ_1=ZZ_1+NI_BHBH
     &	(TEST_X,TEST_Y,Z_M,TEST_LTANPHI,TEST_B,TEST_ALPHA,EXP_X,EXP_Y,
     &	Z_N,EXP_LTANPHI,EXP_B,EXP_ALPHA,SUBREAC_TYPE,K,D,W,EPSILON,
     &	RHO_MAX)			  	
						CASE ('R')
*						  CASE BHRH
						CASE ('L')
*						  CASE BHLH
					    CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT 
					CASE ('F')
*			          CASES FH_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FHFH
						  IF (ABS(ALPHA_MN)==0. 
     &							  .OR. ABS(ALPHA_MN)==2*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE*
*	       				    LOCATE LOOKUP TABLE
						    HEAD => Z_LOOKUP_HH_0 (SEC_M, SEC_N,REGION_M)
*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH(Z_LOOKUP_HH_0 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N
     &                              , L_M, L_N, B_M, B_N, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND
CX                              FOUND =.FALSE.
						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT

!						     MB_L=-1
!						     MB_U=1
!						     SGN_MB=2
!						     NB_L=-1
!						     NB_U=1
!						     SGN_NB=2
!						     ML_L=-1
!						     ML_U=1
!						     SGN_ML=1
!						     NL_L=-1
!						     NL_U=1
!						     SGN_NL=1
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1)

							 ZZ_1=ZZ_1+FHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						   NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
*							  ADD ENTRY TO LOOKUP TABLE

						HEAD => Z_LOOKUP_HH_0 (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT(Z_LOOKUP_HH_0 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                     L_M, L_N, B_M, B_N,ZZ_1)
C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSEIF (ABS(ALPHA_MN)==PI 
     &							  .OR. ABS(ALPHA_MN)==3*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
					HEAD => Z_LOOKUP_HH_PI (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH(Z_LOOKUP_HH_0
     &				(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND

CX                              FOUND =.FALSE.
						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1)

							 ZZ_1=ZZ_1+FHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
*							  ADD ENTRY TO LOOKUP TABLE

					 HEAD => Z_LOOKUP_HH_0 (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT(Z_LOOKUP_HH_0
     &				(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSE IF (     ALPHA_MN==PI/2 
     &							   .OR. ALPHA_MN==-3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*

				HEAD => Z_LOOKUP_HH_90_PLUS (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH(Z_LOOKUP_HH_90_PLUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)	

CX                           FOUND = .FALSE.
						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
	                    
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1)


							 ZZ_1=ZZ_1+FHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)

*							  ADD ENTRY TO LOOKUP TABLE
				HEAD => Z_LOOKUP_HH_90_PLUS (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT(Z_LOOKUP_HH_90_PLUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N,
     &                                ZZ_1 )

							END IF

*---------------------------------------------
*---------------------------------------------
						  ELSE IF (     ALPHA_MN== -PI/2 
     &							   .OR. ALPHA_MN== 3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*

				  HEAD => Z_LOOKUP_HH_90_MINUS (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=( X_M*SIN_ALPHA_MN-Y_N)
	  CALL TABLE_SEARCH(Z_LOOKUP_HH_90_MINUS(SEC_M, SEC_N, REGION_M), 
     &							DDX, DDY, Z_M, Z_N,
     &                            L_M, L_N, B_M, B_N, ZE, FOUND)	
CX	                   FOUND = .FALSE.

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1)


							 ZZ_1=ZZ_1+FHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				 NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &					NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				   STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				   STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)

*							  ADD ENTRY TO LOOKUP TABLE
				HEAD => Z_LOOKUP_HH_90_MINUS (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT(Z_LOOKUP_HH_90_MINUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N,
     &                                ZZ_1)

							END IF

*---------------------------------------------
						  ELSE
*			Arbitrary Angle	ALPHA_MN <>0,PI/2,PI..., COMPUTE THE REACTION 
*			ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1)

							 ZZ_1=ZZ_1+FHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						   NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF
						CASE ('B')
*						  CASE FHBH
		   				  ZZ_1=ZZ_1+NI_FHBH
     &	(TEST_X,TEST_Y,Z_M,TEST_LTANPHI,TEST_B,TEST_ALPHA,EXP_X,EXP_Y,
     &	Z_N,EXP_LTANPHI,EXP_B,EXP_ALPHA,SUBREAC_TYPE,K,D,W,EPSILON,
     &	RHO_MAX)			  	
						CASE ('L')
*						  CASE FHLH
* USING LHFH ROUTINE DUE TO THE RECIPROCITY
						  IF (ABS(ALPHA_MN)==0.) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*							DUE TO SYMMETRY, CAN USE RHFH TABLE WITH DX REPLACED BY -DX INSTEAD
*

*	       				    LOCATE LOOKUP TABLE
                         HEAD => Z_LOOKUP_RFHH(SEC_M, SEC_N ,REGION_M)
                              
*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=-(X_N-X_M)
							DDY=ABS(Y_N-Y_M)
						    CALL TABLE_SEARCH(Z_LOOKUP_RFHH
     &				(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=-(X_N-X_M)
							DDY=ABS(Y_N-Y_M)
						    CALL TABLE_SEARCH(Z_LOOKUP_RFHH
     &					(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)

							  ZZ_TEMP=YHFH(X_N,Y_N,Z_N,
     &						    -L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,-ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,1,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							-SIN_ALPHA_MN,
     &				 NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						  NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &				   STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &				   STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)

*							  ADD ENTRY TO LOOKUP TABLE
                           HEAD => Z_LOOKUP_RFHH(SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT(Z_LOOKUP_RFHH
     &    				(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M,
     &                              Z_N, L_M, L_N, B_M, B_N, ZZ_TEMP)

							  ZZ_1=ZZ_1+ZZ_TEMP
							END IF
						  ELSE
*						ALPHA_MN NOT 0,	COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)

							  ZZ_1=ZZ_1+YHFH(X_N,Y_N,Z_N,
     &						    -L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,-ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,1,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							-SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						   NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF

						CASE ('R')
*						  CASE FHRH
* USING RHFH ROUTINE DUE TO RECIPROCITY
						  IF (ABS(ALPHA_MN)==0) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
                        HEAD => Z_LOOKUP_RFHH (SEC_M, SEC_N, REGION_M)
*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=X_N-X_M
							DDY=ABS(Y_N-Y_M)
						    CALL TABLE_SEARCH(Z_LOOKUP_RFHH 
     &				    (SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)


						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)

							  ZZ_TEMP=YHFH(X_N,Y_N,Z_N,
     &						    L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,-ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							-SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						  NUM_GREEN_HH(LAYER_N,LAYER_M, REGION_M),
     &				   STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M, REGION_M),
     &				   STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N, REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)

*							  ADD ENTRY TO LOOKUP TABLE
                           HEAD => Z_LOOKUP_RFHH(SEC_M, SEC_N, REGION_M)

						      CALL TABLE_INSERT(Z_LOOKUP_RFHH
     &					 (SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M,
     &                             Z_N, L_M, L_N, B_M, B_N, ZZ_TEMP)

							  ZZ_1=ZZ_1+ZZ_TEMP
							END IF
						  ELSE
*		ALPHA_MN NOT 0,	COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)

							  ZZ_1=ZZ_1+YHFH(X_N,Y_N,Z_N,
     &						    L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,-ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							-SIN_ALPHA_MN,
     &				NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &					     NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &				  STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &				  STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('L')
*			          CASES LH_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE LHFH
						  IF (ABS(ALPHA_MN)==0.) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*							DUE TO SYMMETRY, CAN USE RHFH TABLE WITH DX REPLACED BY -DX INSTEAD
**	       				    LOCATE LOOKUP TABLE
                         HEAD => Z_LOOKUP_RFHH (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=-(X_M-X_N)
							DDY=ABS(Y_M-Y_N)
						    CALL TABLE_SEARCH(Z_LOOKUP_RFHH 
     &				   (SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)
	
							  ZZ_TEMP=YHFH(X_M,Y_M,Z_M,
     &						    -L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,1,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
*							  ADD ENTRY TO LOOKUP TABLE

                          HEAD => Z_LOOKUP_RFHH (SEC_M, SEC_N,REGION_M)

						      CALL TABLE_INSERT(Z_LOOKUP_RFHH 
     &					     (SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M,
     &                               Z_N, L_M, L_N, B_M, B_N, ZZ_TEMP)

							  ZZ_1=ZZ_1+ZZ_TEMP
							END IF
						  ELSE
*						ALPHA_MN NOT 0,	COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)
	
							  ZZ_1=ZZ_1+YHFH(X_M,Y_M,Z_M,
     &						    -L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,1,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF

					    CASE ('L')
*						  CASE LHLH
						  IF (ABS(ALPHA_MN)==0.) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*							BECAUSE OF THE SYMMETRY, CAN USE THE RHRH TABLE
*
*	       				    LOCATE LOOKUP TABLE
                          HEAD => Z_LOOKUP_RRHH (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=ABS(X_M-X_N)
							DDY=ABS(Y_M-Y_N)
						    CALL TABLE_SEARCH(Z_LOOKUP_RRHH 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                             L_M, L_N, B_M, B_N, ZE, FOUND)

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*	ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

							  ZZ_TEMP=YHYH(X_M,Y_M,Z_M,
     &						    -L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				 NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)

*							  ADD ENTRY TO LOOKUP TABLE

                          HEAD => Z_LOOKUP_RRHH (SEC_M, SEC_N,REGION_M)

						      CALL TABLE_INSERT(Z_LOOKUP_RRHH 
     &						(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M,
     &                              Z_N, L_M, L_N, B_M, B_N, ZZ_TEMP)

							  ZZ_1=ZZ_1+ZZ_TEMP
							END IF
						  ELSE
*		ALPHA_MN NOT 0,	COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

							  ZZ_1=ZZ_1+YHYH(X_M,Y_M,Z_M,
     &						    -L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &					   	NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				    STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				    STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF
					    CASE ('R')
*						  CASE LHRH
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

						  ZZ_1=ZZ_1+YHYH(X_M,Y_M,Z_M,
     &						    -L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,1,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						CASE ('B')
*						  CASE LHBH
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('R')
*			          CASES RH_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE RHFH
						  IF (ABS(ALPHA_MN)==0) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*

*	       				    LOCATE LOOKUP TABLE
                        HEAD => Z_LOOKUP_RFHH (SEC_M, SEC_N, REGION_M)


*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=X_M-X_N
							DDY=ABS(Y_M-Y_N)
						    CALL TABLE_SEARCH(Z_LOOKUP_RFHH 
     &				(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, B_M, B_N, ZE, FOUND)

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)
	
							  ZZ_TEMP=YHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
*							  ADD ENTRY TO LOOKUP TABLE

                         HEAD => Z_LOOKUP_RFHH (SEC_M, SEC_N, REGION_M)

						      CALL TABLE_INSERT(Z_LOOKUP_RFHH 
     &					     (SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M,
     &                             Z_N, L_M, L_N, B_M, B_N, ZZ_TEMP)

							  ZZ_1=ZZ_1+ZZ_TEMP
							END IF
						  ELSE
*						ALPHA_MN NOT 0,	COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)
	
							  ZZ_1=ZZ_1+YHFH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				 NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						  NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				   STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				   STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF

					    CASE ('L')
*						  CASE RHLH
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

						  ZZ_1=ZZ_1+YHYH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,1,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
					    CASE ('R')
*						  CASE RHRH
						  IF (ABS(ALPHA_MN)==0.) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
                         HEAD => Z_LOOKUP_RRHH (SEC_M, SEC_N, REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=ABS(X_M-X_N)
							DDY=ABS(Y_M-Y_N)
						    CALL TABLE_SEARCH(Z_LOOKUP_RRHH 
     &				(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                             L_M, L_N, B_M, B_N, ZE, FOUND)

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
	
							  ZZ_TEMP=YHYH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &				  NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &						NUM_GREEN_HH(LAYER_M,LAYER_N, REGION_M),
     &				   STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N, REGION_M),
     &				   STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N, REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)

*							  ADD ENTRY TO LOOKUP TABLE

                          HEAD => Z_LOOKUP_RRHH (SEC_M, SEC_N,REGION_M)

						      CALL TABLE_INSERT(Z_LOOKUP_RRHH 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M,
     &                              Z_N, L_M, L_N, B_M, B_N, ZZ_TEMP)

							  ZZ_1=ZZ_1+ZZ_TEMP
							END IF
						  ELSE
*						ALPHA_MN NOT 0,	COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
	
							  ZZ_1=ZZ_1+YHYH(X_M,Y_M,Z_M,
     &						    L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,ALPHA_MN,
     &						    K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,
     &							NB_U,SGN_NB, ML_L,ML_U,SGN_ML,
     &							NL_L,NL_U,SGN_NL,0,RHO_MAX,RHO_CENTER,
     &							THETA_MN_CENTER,COS_ALPHA_MN,
     &							SIN_ALPHA_MN,
     &							NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),
     &							NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &						    STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &						    STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &							MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &							STORE_KRHO_TE(0,REGION_M),GND_LOSSY)
						  END IF

						CASE ('B')
*						  CASE RHBH
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
			    CASE ('V')
*			      CASES _H_V
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*			          CASES FH_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FHFV

	                    T_L=-1
	                    T_U=1
						  E_L=-1
						  E_U=1
						  ZZ_1=ZZ_1+YVFH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
						CASE ('L')
*						  CASE FHLV
                          T_L=-1
	                    T_U=1
						  E_L=0
						  E_U=-1

						  ZZ_1=ZZ_1+YVFH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
						CASE ('R')
*						  CASE FHRV
                          T_L=-1
						  T_U=1   
						  E_L=0
						  E_U=1

						  ZZ_1=ZZ_1+YVFH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('L')
*			          CASES LH_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE LHFV

						  T_L=0
						  T_U=-1
						  E_L=-1
						  E_U=1
						  ZZ_1=ZZ_1+YVYH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
          			    CASE ('L')
*						  CASE LHLV
						  T_L=0
						  T_U=-1
						  E_L=0
						  E_U=-1

						  ZZ_1=ZZ_1+YVYH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
					    CASE ('R')
*						  CASE LHRV
						  T_L=0
						  T_U=-1
						  E_L=0
						  E_U=1

						  ZZ_1=ZZ_1+YVYH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('R')
*			          CASES RH_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE RHFV

						  T_L=0
						  T_U=1
						  E_L=-1
						  E_U=1
						  
						  ZZ_1=ZZ_1+YVYH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
					    CASE ('L')
*						  CASE RHLV
						  T_L=0
						  T_U=1
						  E_L=0
						  E_U=-1

						  ZZ_1=ZZ_1+YVYH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
					    CASE ('R')
*						  CASE RHRV
						  T_L=0
						  T_U=1
						  E_L=0
						  E_U=1

						  ZZ_1=ZZ_1+YVYH(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,T_L,T_U,
     &                        RHO_MAX,CURR_FLAG_N,ALPHA_MN_EQ0)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
*----------------------------------ADD BY WANG K TYPE CELL JUDGE BEGINS -------------------
			    CASE ('K')
*			      CASES _H_K
			      SELECT CASE (REAC_TYPE(III)(1:1))
	                CASE ('L')
*					  CASES LH_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE LHFK
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

	              L_N2=CELL_Z(J)%L_OR_TANPHI(2)

	                   ZZ_1 = ZZ_1 +
     &            FKYH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M,-L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &  -SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),		
     &					 NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2)


						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT


	                CASE ('R')
*					  CASES RH_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE RHFK

				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

	                          L_N2=CELL_Z(J)%L_OR_TANPHI(2)

	                   ZZ_1 = ZZ_1 +
     &            FKYH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M, L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	-SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),			
     &					  NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2)



						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT


					CASE ('F')
*					  CASES FH_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FHFK

	                          L_N2=CELL_Z(J)%L_OR_TANPHI(2)

						  IF (ABS(ALPHA_MN)==0. 
     &							  .OR. ABS(ALPHA_MN)== 2*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
						HEAD1 => Z_LOOKUP_KH_0 (SEC_N, SEC_M,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_0
     &				 (SEC_N, SEC_M,REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                             L_N, L_M, L_N2,
     &                             B_N, B_M, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND

CX                              FOUND =.FALSE.
						    IF (FOUND) THEN
        
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*  ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	-SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),			
     &						NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2)
*							  ADD ENTRY TO LOOKUP TABLE

					HEAD1 => Z_LOOKUP_KH_0 (SEC_N, SEC_M, REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_0
     &				(SEC_N, SEC_M,REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                              L_N, L_M, L_N2,  B_N, B_M,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSE IF (ABS(ALPHA_MN)==PI 
     &							  .OR. ABS(ALPHA_MN)== 3*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
				 HEAD1 => Z_LOOKUP_KH_PI (SEC_N, SEC_M, REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_PI
     &				(SEC_N, SEC_M, REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                             L_N, L_M, L_N2,
     &                             B_N, B_M, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND

CX                              FOUND =.FALSE.
						    IF (FOUND) THEN
        
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &  -SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),		
     &						NUM_GREEN_HH(LAYER_N,LAYER_M, REGION_M),
     &				 STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M, REGION_M),
     &				 STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M, REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2)
*							  ADD ENTRY TO LOOKUP TABLE

					HEAD1 => Z_LOOKUP_KH_PI (SEC_N, SEC_M, REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_PI
     &				(SEC_N, SEC_M, REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                              L_N, L_M, L_N2,  B_N, B_M,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSE IF (     ALPHA_MN==PI/2 
     &							   .OR. ALPHA_MN==-3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*                                FILL THE 90-MINUS TABLE    
 				 HEAD1 => Z_LOOKUP_KH_90_MINUS (SEC_N, SEC_M, REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_90_MINUS 
     &				  (SEC_N, SEC_M, REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                              L_N, L_M, L_N2, 
     &                             B_N, B_M, ZE, FOUND)	

CX                           FOUND = .FALSE.

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
	                    
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	-SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),			
     &						NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2 )

*							  ADD ENTRY TO LOOKUP TABLE
				 HEAD1 => Z_LOOKUP_KH_90_MINUS (SEC_N, SEC_M,REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_90_MINUS 
     &					(SEC_N, SEC_M,REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                             L_N, L_M, L_N2, B_N, B_M,
     &                                ZZ_1 )

							END IF

*---------------------------------------------
*---------------------------------------------
						  ELSE IF (     ALPHA_MN== -PI/2 
     &							   .OR. ALPHA_MN== 3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*                                     FILL THE 90_PLUS TABLE
				HEAD1 => Z_LOOKUP_KH_90_PLUS (SEC_N, SEC_M, REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_90_PLUS 
     &					(SEC_N, SEC_M, REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                              L_N, L_M, L_N2,
     &                             B_N, B_M, ZE, FOUND)	
CX	                   FOUND = .FALSE.


						    IF (FOUND) THEN
C      PRINT *, FOUND
C	PAUSE
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	-SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2 )

*							  ADD ENTRY TO LOOKUP TABLE
				HEAD1 => Z_LOOKUP_KH_90_PLUS (SEC_N, SEC_M,REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_90_PLUS 
     &				(SEC_N, SEC_M,REGION_M), -DDX, -DDY, Z_N, Z_M,
     &                              L_N, L_M,L_N2, B_N, B_M,
     &                                ZZ_1)

							END IF

*---------------------------------------------
						  ELSE
*							ALPHA_MN <>0,PI/2,PI..., COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE

				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_N,Y_N,Z_N,-L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,
     &			-ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &   -SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_N,LAYER_M,REGION_M),
     &				   STORE_GREEN_HH_TM(1,LAYER_N,LAYER_M,REGION_M),
     &				   STORE_GREEN_HH_TE(1,LAYER_N,LAYER_M,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_N2 )
	                     
						 END IF

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
*---------------------------------ADDED BY WANG K TYPE CELL JUDGE ENDS------------

				CASE DEFAULT
				  STOP '***REACTION TYPE MISMATCH IN MAIN***'
			  END SELECT
		    CASE ('V')
*			  CASES _V__
			  SELECT CASE (REAC_TYPE(III)(4:4))
			    CASE ('H')
*			      CASES _V_H
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*			          CASES FV_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FVFH

						  T_L=-1
						  T_U=1
	                      E_L=-1
	                      E_U=1
						  ZZ_1=ZZ_1+YVFH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX,CURR_FLAG_M,ALPHA_MN_EQ0)
						CASE ('L')
*						  CASE FVLH

						  T_L=-1
						  T_U=1
						  E_L=0
						  E_U=-1
						  ZZ_1=ZZ_1+YVYH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
						CASE ('R')
*						  CASE FVRH

						  T_L=-1
						  T_U=1
						  E_L=0
						  E_U=1
						  ZZ_1=ZZ_1+YVYH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('L')
*			          CASES LV_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE LVFH
						  T_L=0
						  T_U=-1
	                      E_L=-1
	                      E_U=1

						  ZZ_1=ZZ_1+YVFH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
					    CASE ('L')
*						  CASE LVLH
						  T_L=0
						  T_U=-1
						  E_L=0
						  E_U=-1

						  ZZ_1=ZZ_1+YVYH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
					    CASE ('R')
*						  CASE LVRH
						  T_L=0
						  T_U=-1
						  E_L=0
						  E_U=1

						  ZZ_1=ZZ_1+YVYH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('R')
*			          CASES RV_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE RVFH
						  T_L=0
						  T_U=1
	                      E_L=-1
	                      E_U=1

						  ZZ_1=ZZ_1+YVFH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
					    CASE ('L')
*						  CASE RVLH
						  T_L=0
						  T_U=1
						  E_L=0
						  E_U=-1

						  ZZ_1=ZZ_1+YVYH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
					    CASE ('R')
*						  CASE RVRH
						  T_L=0
						  T_U=1
						  E_L=0
						  E_U=1

						  ZZ_1=ZZ_1+YVYH(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_M, ALPHA_MN_EQ0)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT

			    CASE ('V')
*			      CASES _V_V
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*			          CASES FV_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FVFV
						  T_L=-1
						  T_U=1
						  E_L=-1
						  E_U=1

						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)

						CASE ('L')
*						  CASE FVLV
						  T_L=-1
						  T_U=1
						  E_L=0
						  E_U=-1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)

						CASE ('R')
*						  CASE FVRV
						  T_L=-1
						  T_U=1
						  E_L=0
						  E_U=1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('L')
*			          CASES LV_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE LVFV
						  T_L=0
						  T_U=-1
						  E_L=-1
						  E_U=1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
					    CASE ('L')
*						  CASE LVLV
						  T_L=0
						  T_U=-1
						  E_L=0
						  E_U=-1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
					    CASE ('R')
*						  CASE LVRV
						  T_L=0
						  T_U=-1
						  E_L=0
						  E_U=1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE ('R')
*			          CASES RV_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
					    CASE ('F')
*						  CASE RVFV
						  T_L=0
						  T_U=1
						  E_L=-1
						  E_U=1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
					    CASE ('L')
*						  CASE RVLV
						  T_L=0
						  T_U=1
						  E_L=0
						  E_U=-1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
					    CASE ('R')
*						  CASE RVRV
						  T_L=0
						  T_U=1
						  E_L=0
						  E_U=1
						  ZZ_1=ZZ_1+YVYV(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,E_L,E_U,
     &					   CURR_FLAG_M,CURR_FLAG_N)
						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
*------------------------------ADDED BY WANG K TYPE CELL JUDGE BEGINS---------
			    CASE ('K')
*			      CASES _V_K
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*					  CASES FV_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FVFK

						  T_L=-1
						  T_U= 1
	                      L_N2= CELL_Z(J)%L_OR_TANPHI(2)

						  ZZ_1=ZZ_1+
     &                        YVFK(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,
     &                        RHO_MAX, CURR_FLAG_M, 
     &                        L_N2)

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT

					CASE ('L')
*					  CASES LV_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE LVFK

						  T_L= 0
						  T_U=-1
	                      L_N2= CELL_Z(J)%L_OR_TANPHI(2)

						  ZZ_1=ZZ_1+
     &                        YVFK(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,
     &                        RHO_MAX, CURR_FLAG_M, 
     &                        L_N2)

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT

					CASE ('R')
*					  CASES RV_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE RVFK
						  T_L= 0
						  T_U= 1
	                      L_N2= CELL_Z(J)%L_OR_TANPHI(2)

						  ZZ_1=ZZ_1+
     &                        YVFK(X_M,Y_M,Z_M,ALPHA_MN,
     &						L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,T_L,T_U,
     &                        RHO_MAX, CURR_FLAG_M, 
     &                        L_N2)

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
*------------------------------ADDED BY WANG K TYPE CELL JUDGE ENDS----------
				CASE DEFAULT
				  STOP '***REACTION TYPE MISMATCH IN MAIN***'
			  END SELECT
*---------------------------ADDED BY WANG K TYPE CELL JUDGE BEGINS----------------
              CASE ('K')   
*			  CASES _K__			          
			  SELECT CASE (REAC_TYPE(III)(4:4))
			    CASE ('K')
*			      CASES _K_K
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*					  CASES FK_K
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FKFK

	                          L_M2=CELL_Z(I)%L_OR_TANPHI(2)
	                          L_N2=CELL_Z(J)%L_OR_TANPHI(2)

						  IF (ABS(ALPHA_MN)==0. 
     &							  .OR. ABS(ALPHA_MN)==2*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE*
*	       				    LOCATE LOOKUP TABLE
					HEAD2 => Z_LOOKUP_KK_0 (SEC_M, SEC_N,REGION_M)
*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH_2(Z_LOOKUP_KK_0 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                             L_M, L_N, L_M2, L_N2, 
     &                             B_M, B_N, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND
CX                              FOUND =.FALSE.
						    IF (FOUND) THEN        
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
                 

	                   ZZ_1 = ZZ_1 +
     &               FKFK(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2, L_N2 )
*							  ADD ENTRY TO LOOKUP TABLE

					HEAD2 => Z_LOOKUP_KK_0 (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT_2(Z_LOOKUP_KK_0 
     &				(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, L_N2, B_M, B_N,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSEIF (ABS(ALPHA_MN)==PI 
     &							  .OR. ABS(ALPHA_MN)==3*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
					HEAD2 => Z_LOOKUP_KK_PI (SEC_M, SEC_N, REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH_2(Z_LOOKUP_KK_PI 
     &					(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                             L_M, L_N, L_M2, L_N2, 
     &                             B_M, B_N, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND

CX                              FOUND =.FALSE.
						    IF (FOUND) THEN
        
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
                 

	                   ZZ_1 = ZZ_1 +
     &               FKFK(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &					NUM_GREEN_HH(LAYER_M,LAYER_N, REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N, REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N, REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2, L_N2 )
*							  ADD ENTRY TO LOOKUP TABLE

					HEAD2 => Z_LOOKUP_KK_PI (SEC_M, SEC_N, REGION_M)

						    CALL TABLE_INSERT_2(Z_LOOKUP_KK_PI 
     &				(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, L_N2, B_M, B_N,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSE IF (     ALPHA_MN==PI/2 
     &							   .OR. ALPHA_MN==-3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*

				HEAD2 => Z_LOOKUP_KK_90_PLUS (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=ABS(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=ABS( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH_2(Z_LOOKUP_KK_90_PLUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, L_N2,
     &                             B_M, B_N, ZE, FOUND)	

CX                          FOUND = .FALSE.
						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
	                    
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
                 
	                   ZZ_1 = ZZ_1 +
     &               FKFK(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2, L_N2 )

*							  ADD ENTRY TO LOOKUP TABLE
				HEAD2 => Z_LOOKUP_KK_90_PLUS (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT_2(Z_LOOKUP_KK_90_PLUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, L_N2, B_M, B_N,
     &                                ZZ_1 )

							END IF
*---------------------------------------------
*---------------------------------------------
						  ELSE IF (     ALPHA_MN== -PI/2 
     &							   .OR. ALPHA_MN== 3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
				HEAD2 => Z_LOOKUP_KK_90_MINUS (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=ABS(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=ABS( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH_2(Z_LOOKUP_KK_90_MINUS 
     &				(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, L_N2,
     &                             B_M, B_N, ZE, FOUND)	
CX	                   FOUND = .FALSE.

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
                 
	                   ZZ_1 = ZZ_1 +
     &               FKFK(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2, L_N2 )

*							  ADD ENTRY TO LOOKUP TABLE
				HEAD2 => Z_LOOKUP_KK_90_MINUS (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT_2(Z_LOOKUP_KK_90_MINUS 
     &				(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N,L_M2, L_N2, B_M, B_N,
     &                                ZZ_1)

							END IF

*---------------------------------------------
						  ELSE
*							ALPHA_MN <>0,PI/2,PI..., COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)
                 
	                          L_M2=CELL_Z(I)%L_OR_TANPHI(2)
	                          L_N2=CELL_Z(J)%L_OR_TANPHI(2)


	                   ZZ_1 = ZZ_1 +
     &               FKFK(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &					NUM_GREEN_HH(LAYER_M,LAYER_N, REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N, REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N, REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2, L_N2 )

                           END IF

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT


			    CASE ('H')
*			      CASES _K_H
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*					  CASES FK_H
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FKFH


	                          L_M2=CELL_Z(I)%L_OR_TANPHI(2)

						  IF (ABS(ALPHA_MN)==0. 
     &							  .OR. ABS(ALPHA_MN)== 2*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
					HEAD1 => Z_LOOKUP_KH_0 (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_0
     &				(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                             L_M, L_N, L_M2,
     &                             B_M, B_N, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND

CX                              FOUND =.FALSE.
						    IF (FOUND) THEN
        
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2 )
*							  ADD ENTRY TO LOOKUP TABLE

						HEAD1 => Z_LOOKUP_KH_0 (SEC_M, SEC_N,REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_0
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2,  B_M, B_N,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSE IF (ABS(ALPHA_MN)==PI 
     &							  .OR. ABS(ALPHA_MN)== 3*PI) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*
*	       				    LOCATE LOOKUP TABLE
						    HEAD1 => Z_LOOKUP_KH_PI (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(X_M*COS_ALPHA_MN-X_N)
							DDY=(Y_M*COS_ALPHA_MN-Y_N)

						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_PI
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                            L_M, L_N, L_M2,
     &                             B_M, B_N, ZE, FOUND)	
C                             PRINT *, "FOUND FLAG=", FOUND

CX                             FOUND =.FALSE.
						    IF (FOUND) THEN
        
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2 )
*							  ADD ENTRY TO LOOKUP TABLE

				HEAD1 => Z_LOOKUP_KH_PI (SEC_M, SEC_N, REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_PI
     &					(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                             L_M, L_N, L_M2,  B_M, B_N,
     &                                ZZ_1)

C                             PRINT *, "WRITE TABLE....."

							END IF
*---------------------------------------------
						  ELSE IF (     ALPHA_MN==PI/2 
     &							   .OR. ALPHA_MN==-3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*

				HEAD1 => Z_LOOKUP_KH_90_PLUS (SEC_M, SEC_N,REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_90_PLUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, 
     &                             B_M, B_N, ZE, FOUND)	

CX                           FOUND = .FALSE.

						    IF (FOUND) THEN
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
	                    
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)


	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	 SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),			
     &							NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &					STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2 )

*							  ADD ENTRY TO LOOKUP TABLE
				HEAD1 => Z_LOOKUP_KH_90_PLUS (SEC_M, SEC_N, REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_90_PLUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2, B_M, B_N,
     &                                ZZ_1 )

							END IF

*---------------------------------------------
*---------------------------------------------
						  ELSE IF (     ALPHA_MN== -PI/2 
     &							   .OR. ALPHA_MN== 3*PI/2 ) THEN
*						    CHECK TO SEE IF THIS IMPEDANCE ELEMENT HAS ALREADY
*						    BEEN COMPUTED AND STORED IN THE LOOK-UP TABLE
*

				 HEAD1 => Z_LOOKUP_KH_90_MINUS (SEC_M, SEC_N, REGION_M)

*						    CHECK TO SEE IF ENTRY EXISTS IN LOOKUP TABLE
							DDX=(-Y_M*SIN_ALPHA_MN-X_N)
							DDY=( X_M*SIN_ALPHA_MN-Y_N)
						    CALL TABLE_SEARCH_1(Z_LOOKUP_KH_90_MINUS 
     &					(SEC_M, SEC_N,REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N, L_M2,
     &                             B_M, B_N, ZE, FOUND)	
CX	                   FOUND = .FALSE.


						    IF (FOUND) THEN
C      PRINT *, FOUND
C	PAUSE
*						      ENTRY EXISTS IN LOOKUP TABLE
						      ZZ_1=ZZ_1+ZE
C							  PRINT*, 'USE TABLE'
						    ELSE
*							 ENTRY DOES NOT EXIST IN LOOKUP TABLE, SO CALCULATE IT
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)

	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &   SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2 )

*							  ADD ENTRY TO LOOKUP TABLE
			HEAD1 => Z_LOOKUP_KH_90_MINUS (SEC_M, SEC_N, REGION_M)

						    CALL TABLE_INSERT_1(Z_LOOKUP_KH_90_MINUS 
     &					(SEC_M, SEC_N, REGION_M), DDX, DDY, Z_M, Z_N,
     &                              L_M, L_N,L_M2, B_M, B_N,
     &                                ZZ_1)

							END IF

*---------------------------------------------
						  ELSE
*							ALPHA_MN <>0,PI/2,PI..., COMPUTE THE REACTION ELEMENT WITHOUT A LOOKUP TABLE
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,2)

	                   ZZ_1 = ZZ_1 +
     &               FKFH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &  SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),			
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2 )

	                  END IF
*---------------------------------------------
						CASE ('L')
*						  CASE FKLH
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

	                          L_M2=CELL_Z(I)%L_OR_TANPHI(2)

	                   ZZ_1 = ZZ_1 +
     &            FKYH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,-L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &  SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),			
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2)
 

						CASE ('R')
*						  CASE FKRH
				    ASSIGNCOEFF(MB_L,MB_U,SGN_MB,NB_L,NB_U,SGN_NB,
     &					   ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,3)

	                          L_M2=CELL_Z(I)%L_OR_TANPHI(2)

	                   ZZ_1 = ZZ_1 +
     &            FKYH(X_M,Y_M,Z_M,-L_M,B_M,X_N,Y_N,Z_N,L_N,B_N,
     &			ALPHA_MN,K,D,W,EPSILON,MB_L,MB_U,SGN_MB,NB_L,NB_U,
     &		   SGN_NB,ML_L,ML_U,SGN_ML,NL_L,NL_U,SGN_NL,1,
     &			RHO_MAX,RHO_CENTER,THETA_MN_CENTER,COS_ALPHA_MN,
     &	SIN_ALPHA_MN,NUM_GREEN_ARRAY,NUM_KRHO,STORE_KRHO(0,REGION_M),				
     &						NUM_GREEN_HH(LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TM(1,LAYER_M,LAYER_N,REGION_M),
     &				 STORE_GREEN_HH_TE(1,LAYER_M,LAYER_N,REGION_M),
     &			MU,SIGMAG,STORE_KRHO_TM(0,REGION_M),
     &            STORE_KRHO_TE(0,REGION_M),GND_LOSSY,L_M2) 

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
			    CASE ('V')
*			      CASES _K_V
			      SELECT CASE (REAC_TYPE(III)(1:1))
					CASE ('F')
*					  CASES FK_V
					  SELECT CASE (REAC_TYPE(III)(3:3))
						CASE ('F')
*						  CASE FKFV
						  E_L=-1
						  E_U= 1

	                      L_M2=CELL_Z(I)%L_OR_TANPHI(2)

							
						  ZZ_1=ZZ_1+
     &  					    YVFK(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_N,
     &                        L_M2)


						CASE ('L')
*						  CASE FKLV
						  E_L= 0
						  E_U=-1

	                      L_M2=CELL_Z(I)%L_OR_TANPHI(2)
							
						  ZZ_1=ZZ_1+
     &  					    YVFK(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_N,
     &                        L_M2)

						CASE ('R')
*						  CASE FKRV


						  E_L= 0
						  E_U= 1

	                      L_M2=CELL_Z(I)%L_OR_TANPHI(2)

							
						  ZZ_1=ZZ_1+
     &  					    YVFK(X_N,Y_N,Z_N,-ALPHA_MN,
     &						L_N,B_N,X_M,Y_M,Z_M,L_M,B_M,K,D,W,EPSILON
     &						,STORE_KRHO,NUM_KRHO,E_L,E_U,
     &                        RHO_MAX, CURR_FLAG_N,
     &                        L_M2)

						CASE DEFAULT
						  STOP '***REACTION TYPE MISMATCH IN MAIN***'
					  END SELECT
					CASE DEFAULT
					  STOP '***REACTION TYPE MISMATCH IN MAIN***'
				  END SELECT
				CASE DEFAULT
				  STOP '***REACTION TYPE MISMATCH IN MAIN***'
			  END SELECT

*---------------------------ADDED BY WANG K TYPE CELL JUDGE ENDS------------------------
			CASE DEFAULT
			  STOP '***REACTION TYPE MISMATCH IN MAIN***'
		  END SELECT		
CXX	      END IF
		END IF
**************************************************
***
**************************************************

		ENDDO     !  KN
		ENDDO     !  KM
		
		ZResult(I,J)=ZZ_1
		ZResult(J,I)=ZZ_1
		
        WRITE(60,*) I,J,REAL(ZZ_1), IMAG(ZZ_1)
        WRITE(*,*)     I,J,REAL(ZZ_1), IMAG(ZZ_1)   					    
         ZZ_1=0
*********---------THE LAST ELEMENT ADDED HERE END.	
*      J LOOP (DO WHILE)	 
	ENDDO
*      I LOOP 
	ENDDO
	
	close(60)
	
	NZ=INDEX-1
*	STOP TIMER FOR MATRIX FILL TIME

      TIME1 = CPSEC()
      
	WRITE(150,*) FREQUENCY, TIME1-TIME0
	WRITE(*,*) 'FREQUENCY=',FREQUENCY
	WRITE(*,*) 'MATRIX FILL TIME (SECONDS): ',TIME1-TIME0
	PRINT *, "NZ=", NZ
	PRINT *, "SPARSITY=", NZ/(NUM_CELL**2+0.D0)


	PRINT*, 'NUM_CELL=',NUM_CELL
	PRINT*, 'FF_COUNT=',FF_COUNT
	PRINT*, 'KK_COUNT=',KK_COUNT
	PRINT*, 'KH_COUNT=',KH_COUNT

******************** LOSSY  **********************************************
	IF(LOSSY==1) THEN

C	PRINT *,MU,SIGMAT
C	PAUSE
	SEGS_T=1
	OPEN(3,FILE='ZZ.TXT',STATUS='UNKNOWN')
C	OPEN(4,FILE='DCR.TXT',STATUS='UNKNOWN')
	DO I=1,NUM_CELL
C	DO I=19,19
	  DO J=I,NUM_CELL
C	  DO J=20,20

	ZZ_SURF=(0.D0,0.D0)

*     LOCATE ZZ(I,J) AND GET VALUE OF ZZ(I,J)
      CALL MATRIX_SEARCH(AA_CG(I), J, ZZ_IJ)

	IF(ABS(I-J)==SEGS_W .OR. I==J) THEN
		ZZ_SURF=SURFACE(CELL_Z(I),CELL_Z(J),FREQUENCY,MU,SIGMAT)
		ZZ_IJ=ZZ_IJ+ZZ_SURF
          
		CALL MATRIX_FILL(AA_CG(I), J, ZZ_IJ)

		IF(I/=J) THEN
C		ZZ(J,I)=ZZ(J,I)+ZZ_SURF
          CALL MATRIX_FILL(AA_CG(J), I, ZZ_IJ)
		END IF
		WRITE(3,*) I, J, REAL(ZZ_SURF),IMAG(ZZ_SURF)
	END IF

	  END DO   !!!!! DO J
	  
	  LEN_C=0
	  DO L=1,CELL_Z(I)%NUM_SUBCELL
		LEN_C=LEN_C+CELL_Z(I)%L_OR_TANPHI(L)
	  END DO
	  IF(CELL_Z(I)%FLAG==1) THEN
		LEN_C=LEN_C*2
	  END IF
	  DCRESIS=SEGS_T*LEN_C/SIGMAT/CELL_Z(I)%B(1)/2.D0/0.00002

	END DO   !!!! DO I

	CLOSE (3)

	END IF
********************* END OF LOSSY  *********************************************
***********************************************************************************
	
	IF (ANS_Y_PARAM=='Y' .OR. ANS_Y_PARAM=='Y') THEN
*FIND Y-PARAMETERS AND S-PARAMETERS FOR DEEMBEDDED LINE PROBLEM

		IF (NUM_PORT==1) THEN

			DO I=1,NUM_CELL
					V(I)=(0.D0,0.D0)
			ENDDO
			
			DO I=1,SS
			  DO J=1,SC(I)
				 V(SOURSE(I,J))=0.D0-S_VALUE(I)
			  END DO
			END DO	 
			
		DO I=1,LDS
		  DO J=1,ZC(I)
              CALL MATRIX_SEARCH(AA_CG(LOAD(I,J)), LOAD(I,J), ZZ_IJ)
              ZZ_IJ = ZZ_IJ + Z_VALUE(I)
              CALL MATRIX_FILL(AA_CG(LOAD(I,J)), LOAD(I,J), ZZ_IJ)
		  END DO
		END DO
		
	TOLER = 1D-6 
      T_CG1 = CPSEC()
*****ZCG_MINE(SET_ERR, N, AAA,AAA_DIAG, B, X, R, TMP, P, Q, Z)
              CALL ZCG_MINE(TOLER, NUM_CELL, AA_CG, ZZ_DIAG, V, II,
     &                      R_CG, TMP_CG, P_CG, Q_CG, Z_CG )
      T_CG2 = CPSEC()

		DELX=CELL_Z(PORT11(1))%X(1)-CELL_Z(PORT1(1))%X(1)

C	PRINT *,DELX

		I_PORT1=(0.D0,0.D0)
		DO I=1,PC1
			I_PORT1=I_PORT1+II(PORT1(I))
C     &			*2.D0*CELL_Z(PORT1(I))%B(1)
C			PRINT *,II(PORT1(I))
		END DO
C	PRINT *,'I ON PORT1',I_PORT1

		I_PORT11=(0.D0,0.D0)
		DO I=1,PC11
			I_PORT11=I_PORT11+II(PORT11(I))
C     &			*2.D0*CELL_Z(PORT11(I))%B(1)
C			PRINT *,II(PORT11(I))
		END DO
C	PRINT *,'I ON PORT2',I_PORT11

		A1F= (I_PORT1*EXP( IM*K*DELX/2.D0)-I_PORT11*
     &		 EXP(-IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)
		B1F=-(I_PORT1*EXP(-IM*K*DELX/2.D0)-I_PORT11*
     &		 EXP( IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)
*COMPUTE S-PARAMETERS	
			S_11=B1F/A1F
			WRITE(115,*) FREQUENCY,REAL(S_11),IMAG(S_11)
			WRITE(*,*) ' S11 = ',S_11
*TOUCHSTONE OUTPUT
			WRITE(155,1551) FREQUENCY,REAL(S_11),IMAG(S_11)
1551			FORMAT(F20.6,2E18.6)

***********----------------- NUM_PORT==2 ----------------------
		ELSE IF (NUM_PORT==2) THEN
*FOREWARD MEASUREMENT
			DO I=1,NUM_CELL
					V(I)=(0.D0,0.D0)
			ENDDO

C	PRINT *,SS/2

			DO I=1,SS/2
			  DO J=1,SC(I)
				 V(SOURSE(I,J))=0.D0-S_VALUE(I)
C				 V(SOURSE(I,J))=0.D0-S_VALUE(I)*1.D9/FREQUENCY
			  END DO
			END DO	 
		DO I=1,LDS
		  DO J=1,ZC(I)
C			ZZ(LOAD(I,J),LOAD(I,J))=ZZ(LOAD(I,J),LOAD(I,J))+Z_VALUE(I)
              CALL MATRIX_SEARCH(AA_CG(LOAD(I,J)), LOAD(I,J), ZZ_IJ)
              ZZ_IJ = ZZ_IJ + Z_VALUE(I)
              CALL MATRIX_FILL(AA_CG(LOAD(I,J)), LOAD(I,J), ZZ_IJ)
		  END DO
		END DO

C		CALL DLSACG(NUM_CELL,ZZ,MAX_CELLS,V,1,II)

      TOLER = 1D-6
      T_CG1 = CPSEC()
              CALL ZCG_MINE(TOLER, NUM_CELL, AA_CG, ZZ_DIAG, V, II,
     &                      R_CG, TMP_CG, P_CG, Q_CG, Z_CG )	
      T_CG2 = CPSEC()

		DELX=CELL_Z(PORT11(1))%X(1)-CELL_Z(PORT1(1))%X(1)

		I_PORT1=(0.D0,0.D0)
		DO I=1,PC1
			I_PORT1=I_PORT1+II(PORT1(I))
C     &			*2.D0*CELL_Z(PORT1(I))%B(1)
C			PRINT *,II(PORT1(I))
		END DO
C	PRINT *,'I ON PORT1',I_PORT1

		I_PORT11=(0.D0,0.D0)
		DO I=1,PC11
			I_PORT11=I_PORT11+II(PORT11(I))
		END DO
C	PRINT *,'I ON PORT2',I_PORT11

		A1F= (I_PORT1*EXP( IM*K*DELX/2.D0)-I_PORT11*
     &		 EXP(-IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)
		B1F=-(I_PORT1*EXP(-IM*K*DELX/2.D0)-I_PORT11*
     &		 EXP( IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)

		I_PORT2=(0.D0,0.D0)
		DO I=1,PC2
			I_PORT2=I_PORT2+II(PORT2(I))
C     &			*2.D0*CELL_Z(PORT2(I))%B(1)
		END DO
C	PRINT *,'I ON PORT2',I_PORT2

		I_PORT22=(0.D0,0.D0)
		DO I=1,PC22
			I_PORT22=I_PORT22+II(PORT22(I))
C     &			*2.D0*CELL_Z(PORT22(I))%B(1)
		END DO

		B2F= (I_PORT2*EXP( IM*K*DELX/2.D0)-I_PORT22*
     &		 EXP(-IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)
		A2F=-(I_PORT2*EXP(-IM*K*DELX/2.D0)-I_PORT22*
     &		 EXP( IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)

*REVERSE MEASUREMENT

			DO I=1,NUM_CELL
					V(I)=(0.D0,0.D0)
			ENDDO
			
			DO I=SS/2+1,SS
			  DO J=1,SC(I)
				 V(SOURSE(I,J))=0.D0-S_VALUE(I)
C				 V(SOURSE(I,J))=0.D0-S_VALUE(I)*1.D9/FREQUENCY
			  END DO
			END DO	 	
C		CALL DLSACG(NUM_CELL,ZZ,MAX_CELLS,V,1,II)
      TOLER = 1D-6
      T_CG3 = CPSEC()
              CALL ZCG_MINE(TOLER, NUM_CELL, AA_CG, ZZ_DIAG, V, II,
     &                      R_CG, TMP_CG, P_CG, Q_CG, Z_CG )	
      T_CG4 = CPSEC()
        PRINT *,  'MATRIX SOLVE TIME BY  MY OWN
     & CONJUGATE GRADIENT (SECONDS): ', T_CG4-T_CG3+T_CG2-T_CG1

C	PRINT *,'REVERSE MATRIX SOLVED'
		I_PORT1=(0.D0,0.D0)
		DO I=1,PC1
			I_PORT1=I_PORT1+II(PORT1(I))
C     &			*2.D0*CELL_Z(PORT1(I))%B(1)
		END DO
C	PRINT *,'I ON PORT1',I_PORT1

		I_PORT11=(0.D0,0.D0)
		DO I=1,PC11
			I_PORT11=I_PORT11+II(PORT11(I))
C     &			*2.D0*CELL_Z(PORT11(I))%B(1)
		END DO
C	PRINT *,'I ON PORT11',I_PORT11

		A1R= (I_PORT1*EXP( IM*K*DELX/2.D0)-I_PORT11*
     &		 EXP(-IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)
		B1R=-(I_PORT1*EXP(-IM*K*DELX/2.D0)-I_PORT11*
     &		 EXP( IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)

		I_PORT2=(0.D0,0.D0)
		DO I=1,PC2
			I_PORT2=I_PORT2+II(PORT2(I))
C     &			*2.D0*CELL_Z(PORT2(I))%B(1)
C			PRINT *,II(PORT2(I))
		END DO
C	PRINT *,'I ON PORT2',I_PORT2

		I_PORT22=(0.D0,0.D0)
		DO I=1,PC22
			I_PORT22=I_PORT22+II(PORT22(I))
C     &			*2.D0*CELL_Z(PORT22(I))%B(1)
C			PRINT *,II(PORT22(I))
		END DO
C	PRINT *,'I ON PORT22',I_PORT22

		B2R= (I_PORT2*EXP( IM*K*DELX/2.D0)-I_PORT22*
     &		 EXP(-IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)
		A2R=-(I_PORT2*EXP(-IM*K*DELX/2.D0)-I_PORT22*
     &		 EXP( IM*K*DELX/2.D0))/2.D0/IM/SIN(K*DELX)

*******************************************************************************

*FILL IN THE MATRIX FOR SOLVING FOR S-PARAMETER		
			DO G=1,4
	 			DO GG=1,4
					A_MATRIX(G,GG)=(0.D0,0.D0)
				END DO
			END DO
		A_MATRIX(1,1)= A1F
		A_MATRIX(1,2)=-A2F
		A_MATRIX(2,1)= A1R
		A_MATRIX(2,2)=-A2R
		A_MATRIX(3,3)= A1F
		A_MATRIX(3,4)=-A2F
		A_MATRIX(4,3)= A1R
		A_MATRIX(4,4)=-A2R
		B_ARRAY(1)=-B1F
		B_ARRAY(2)=-B1R
		B_ARRAY(3)= B2F
		B_ARRAY(4)= B2R

*SOLVE FOR S-PARAMETER
			CALL DLSACG(4,A_MATRIX,16,B_ARRAY,1,S_ARRAY)
			S_11=S_ARRAY(1)
			S_12=S_ARRAY(2)
			S_21=S_ARRAY(3)
			S_22=S_ARRAY(4)
			TIME1=CPSEC()
			WRITE(115,*) FREQUENCY,REAL(S_11),IMAG(S_11)
			WRITE(*,*) ' S11 = ',S_11
			WRITE(116,*) FREQUENCY,REAL(S_21),IMAG(S_21)
			WRITE(*,*) ' S21 = ',S_21
			WRITE(117,*) FREQUENCY,REAL(S_12),IMAG(S_12)
			WRITE(*,*) ' S12 = ',S_12
			WRITE(118,*) FREQUENCY,REAL(S_22),IMAG(S_22)
			WRITE(*,*) ' S22 = ',S_22
			WRITE(*,*) ' '
			WRITE(151,*) FREQUENCY, TIME1-TIME0
*    TOUCHSTONE OUTPUT
		WRITE(155,1552) FREQUENCY,REAL(S_11),IMAG(S_11),REAL(S_21),
     &		IMAG(S_21),REAL(S_12),IMAG(S_12),REAL(S_22),IMAG(S_22)
1552			FORMAT(F20.6,8E18.6)

		END IF
	ENDIF

		OPEN(UNIT=30,FILE="OUTPUT_CURR.DAT",STATUS="REPLACE")
		DO J=1,NUM_CELL
			WRITE(30,*) REAL(II(J)),IMAG(II(J)),ABS(II(J))
		ENDDO
		CLOSE (30)

C	ENDIF
***********************************************************************************
***********************************************************************************
	IF (ANSWER_SWEEP=='Y' .OR. ANSWER_SWEEP=='Y') THEN
		FREQUENCY=FREQUENCY+(FREQ_END-FREQ_START)/FREQ_STEP
	ENDIF
	
      DEALLOCATE (AA_CG)
	DEALLOCATE (R_CG)
	DEALLOCATE (TMP_CG)
	DEALLOCATE (P_CG)
	DEALLOCATE (Q_CG)
	DEALLOCATE  (Z_CG)

	ENDDO



      TIME3 = CPSEC()
      PRINT *, "ALL PROCESS SOLVE TIME IS:" , TIME3-TIME2

	IF (ANSWER_SWEEP=='Y' .OR. ANSWER_SWEEP=='Y') CLOSE (100)

	IF (ANS_Y_PARAM=='Y' .OR. ANS_Y_PARAM=='Y') THEN

		CLOSE (115)
		CLOSE (116)
		CLOSE (117)
		CLOSE (118)

		CLOSE (155)

	ENDIF


		CLOSE (150)
		CLOSE (151)

	IF (ANSWER_GENERATE=='Y' .OR. ANSWER_GENERATE=='Y') THEN

	REPEAT = 1;

	DO WHILE ( REPEAT == 1)
		WRITE(*,*) ' '
		WRITE(*,*) 'ON WHICH LAYER WOULD YOU LIKE TO DISPLAY
     & THE CURRENT:(1,2..)'
		READ(*,*) LAYER_PRINT

		ANSWER_LAYER = ICHAR(LAYER_PRINT) - 48

		WRITE(*,*) 'ON WHICH PLANE WOULD YOU LIKE TO DISPLAY
     & THE CURRENT:(T=TOP/B=BOTTOM)'
		READ(*,*) ANSWER_PLANE


	JSXT_REAL='JSXT_REAL'//LAYER_PRINT//ANSWER_PLANE//'.DAT'
	JSXT_IMAG='JSXT_IMAG'//LAYER_PRINT//ANSWER_PLANE//'.DAT'
	JSYT_REAL='JSYT_REAL'//LAYER_PRINT//ANSWER_PLANE//'.DAT'
	JSYT_IMAG='JSYT_IMAG'//LAYER_PRINT//ANSWER_PLANE//'.DAT'


	OPEN(UNIT=30,FILE=JSXT_REAL)
	OPEN(UNIT=31,FILE=JSXT_IMAG)
	OPEN(UNIT=32,FILE=JSYT_REAL)
	OPEN(UNIT=33,FILE=JSYT_IMAG)

	MICRON = 1.D-6

		IF (ANSWER_PLANE=='T' .OR. ANSWER_PLANE=='T')THEN
	  Z_M = GLOBAL_HEIGHT(ANSWER_LAYER)- GLOBAL_HEIGHT(ANSWER_LAYER-1)

		ELSE
			Z_M = 0 * MICRON
		END IF

	D = GLOBAL_HEIGHT(ANSWER_LAYER)- GLOBAL_HEIGHT(ANSWER_LAYER-1)

      UP_INDICATOR = 1

	IF(1) THEN
	WRITE(*,*) 'WOULD YOU LIKE TO SET THE GROUND CURRENT GRID: (Y/N)?'
		READ(*,*) RESOLUTION

		IF (RESOLUTION=='Y' .OR. RESOLUTION=='Y') THEN

			WRITE(*,*) 'ENTER THE X START POSITION: (IN MICRONS)'
			READ(*,*) X_START
			X_START = X_START * MICRON

			WRITE(*,*) 'ENTER THE X END POSITION: (IN MICRONS)'
			READ(*,*) X_END 
			X_END = X_END * MICRON

			WRITE(*,*) 'ENTER THE NUMBER OF STEPS IN THE X DIRECTION:'
			READ(*,*) X_STEP
			X_STEP = (X_END-X_START)/X_STEP

			WRITE(*,*) 'ENTER THE Y START POSITION: (IN MICRONS)'
			READ(*,*) Y_START
			Y_START = Y_START * MICRON

			WRITE(*,*) 'ENTER THE Y END POSITION: (IN MICRONS)'
			READ(*,*) Y_END
			Y_END = Y_END * MICRON

			WRITE(*,*) 'ENTER THE NUMBER OF STEPS IN THE Y DIRECTION:'
			READ(*,*) Y_STEP
			Y_STEP = (Y_END-Y_START)/Y_STEP

		ELSE
			X_START = X_MIN - L_M - 5*D
			X_END   = X_MAX + L_M + 5*D
			X_STEP  = ((X_MAX + L_M)-(X_MIN - L_M))/60

			Y_START = Y_MIN - B_M - 5*D
			Y_END   = Y_MAX + B_M + 5*D
			Y_STEP  = ((Y_MAX + B_M)-(Y_MIN - B_M))/60

		END IF

	END IF

		WRITE(*,*) 'NOW GENERATING THE GROUND CURRENT DISPLAY '

      DO X_M = X_START, X_END, X_STEP
	DO Y_M = Y_START, Y_END, Y_STEP

       JSXT_TOTAL = 0
	 JSYT_TOTAL = 0 


      DO J= 1, NUM_CELL
		
		SEC_N=CELL_Z(J)%SECTION
		LAYER_N=CELL_Z(J)%LAYER
		FLAG_N=CELL_Z(J)%FLAG
		SUBCELL_N=CELL_Z(J)%NUM_SUBCELL		

	    ALPHA_M = 0

          DO III = 1, SUBCELL_N

	IF(REGION(J,III)==ANSWER_LAYER) THEN

					Z_N = CELL_Z(J)%Z - GLOBAL_HEIGHT(ANSWER_LAYER-1)

					X_N=CELL_Z(J)%X(III)
					Y_N=CELL_Z(J)%Y(III)
					B_N=CELL_Z(J)%B(III)
					L_N=CELL_Z(J)%L_OR_TANPHI(III)
					ALPHA_N=CELL_Z(J)%ALPHA(III)
					CURR_FLAG_N=CELL_Z(J)%CURR_FLAG(III)

				TANDEL = GLOBAL_TANDEL(ANSWER_LAYER)
				EPSILONR = GLOBAL_PERMITTIVITY(ANSWER_LAYER)
				EPSILON = EPSILONO*EPSILONR*(1.D0-IM*TANDEL)
				K = GLOBAL_K(ANSWER_LAYER)


				RHO_MAX=ABS(LOG(RELERR_RHOMN_HH)/DBLE(SQRT((PI/D)**2-K*K)))

							ALPHA_MN= ALPHA_M - ALPHA_N
							IF (ALPHA_MN .EQ. 0.D0) THEN
								ALPHA_MN_EQ0=.TRUE.
							ELSE
								ALPHA_MN_EQ0=.FALSE.
							END IF


	           COS_ALPHA_N = COS(ALPHA_N)
	           SIN_ALPHA_N = SIN(ALPHA_N)
	           IF (ABS(SIN_ALPHA_N) == 1)  COS_ALPHA_N=0.D0
		       IF (ABS(COS_ALPHA_N) == 1)  SIN_ALPHA_N=0.D0  

			REAC_TYPE(III)(1:1)=CELL_Z(J)%LR_CODE(III:III)
			REAC_TYPE(III)(2:2)=CELL_Z(J)%CODE(III:III)

	
		  SELECT CASE (REAC_TYPE(III)(2:2))
		    CASE ('H')
*			  CASES _H
			  SELECT CASE (REAC_TYPE(III)(1:1))
			    CASE ('F')
*			    CASES FH
      
			   E_L=-1
			   E_U=1

     	           JSXE = JSXE_FH(X_M,Y_M,Z_M,ALPHA_MN,X_N,Y_N,
     &				Z_N,L_N,B_N,K,D,W,EPSILON,STORE_KRHO,
     &				NUM_KRHO,E_L,E_U,RHO_MAX)
	           JSYE = 0
	
			    CASE ('R')
*			    CASES RH
      
			   E_L=0
			   E_U=1
     	
	           JSXE = JSXE_YH(X_M,Y_M,Z_M,ALPHA_MN,X_N,Y_N,
     &				Z_N,L_N,B_N,K,D,W,EPSILON,STORE_KRHO,
     &				NUM_KRHO,E_L,E_U,RHO_MAX)
	           JSYE = 0

			    CASE ('L')
*			    CASES LH
      
			   E_L=-1
			   E_U=0
     	
	           JSXE = JSXE_YH(X_M,Y_M,Z_M,ALPHA_MN,X_N,Y_N,
     &				Z_N,L_N,B_N,K,D,W,EPSILON,STORE_KRHO,
     &				NUM_KRHO,E_L,E_U,RHO_MAX)
	           JSYE = 0

				CASE DEFAULT
				  STOP '***REACTION TYPE MISMATCH IN MAIN***'
			  END SELECT


		    CASE ('V')
			IF(1) THEN
*			  CASES _V
			  SELECT CASE (REAC_TYPE(III)(1:1))
			    CASE ('F')
*			    CASES FV      
			   E_L=-1
			   E_U=1
			    CASE ('R')
*			    CASES FV
			   E_L=0
			   E_U=1
			    CASE ('L')
*			    CASES FV
			   E_L=-1
			   E_U=0

				CASE DEFAULT
				  STOP '***REACTION TYPE MISMATCH IN MAIN***'
			  END SELECT

	           JSXE = JSXE_YV(X_M,Y_M,Z_M,ALPHA_MN,X_N,Y_N,
     &				Z_N,L_N,B_N,K,D,W,EPSILON,STORE_KRHO,
     &				NUM_KRHO,E_L,E_U,RHO_MAX, 1)

	           JSYE = JSYE_YV(X_M,Y_M,Z_M,ALPHA_MN,X_N,Y_N,
     &				Z_N,L_N,B_N,K,D,W,EPSILON,STORE_KRHO,
     &				NUM_KRHO,E_L,E_U,RHO_MAX, 1)

			END IF

			CASE DEFAULT
			  STOP '***REACTION TYPE MISMATCH IN MAIN***'
		  END SELECT

                 JSXT = JSXE * COS_ALPHA_N  - JSYE * SIN_ALPHA_N
	           JSYT = JSXE * SIN_ALPHA_N  + JSYE * COS_ALPHA_N

	           JSXT_TOTAL = JSXT_TOTAL + JSXT*II(J)
	           JSYT_TOTAL = JSYT_TOTAL + JSYT*II(J)

	END IF
*         END OF III LOOP
		END DO
C     	END IF
*         END OF J LOOP
       	END DO
	 
	   WRITE(30,*) X_M, Y_M, REAL(JSXT_TOTAL)
		 WRITE(31,*) X_M, Y_M, IMAG(JSXT_TOTAL)
		 WRITE(32,*) X_M, Y_M, REAL(JSYT_TOTAL)
		 WRITE(33,*) X_M, Y_M, IMAG(JSYT_TOTAL)

	END DO
	END DO

	CLOSE(30)
	CLOSE(31)
	CLOSE(32)
	CLOSE(33)
*------------------------------------------------
      IF (1) THEN
                             
		OPEN(UNIT=34,FILE="JS_CONDUCTOR_REAL.DAT",STATUS="REPLACE")
	    OPEN(UNIT=35,FILE="JS_CONDUCTOR_IMAG.DAT",STATUS="REPLACE")
		DO J=1,NUM_CELL
			X_N=CELL_Z(J)%X(1)
			Y_N=CELL_Z(J)%Y(1)
	        B_N=CELL_Z(J)%B(1)
              ALPHA_N=CELL_Z(J)%ALPHA(1)
          
	    COS_ALPHA_N = COS(ALPHA_N)
	    SIN_ALPHA_N = SIN(ALPHA_N)
	    IF (ABS(SIN_ALPHA_N) == 1)  COS_ALPHA_N=0.D0
		IF (ABS(COS_ALPHA_N) == 1)  SIN_ALPHA_N=0.D0  
            
		   DO I= -1, 1, 2
* CHANGE TO GLOBAL COORDINATES
	X_NO = X_N*COS_ALPHA_N - (Y_N+I*B_N)*SIN_ALPHA_N
      Y_NO = X_N*SIN_ALPHA_N + (Y_N+I*B_N)*COS_ALPHA_N


		WRITE(34,*) ' ',X_NO, Y_NO, REAL(II(J))
          WRITE(35,*) ' ',X_NO, Y_NO, IMAG(II(J))
	       ENDDO
		ENDDO
		
		CLOSE (34)
	    CLOSE (35)

      END IF

	WRITE(*,*) ' '
	WRITE(*,*) 'DATA OUTPUT TO THESE FILES:'
	WRITE(*,*) ' '
	WRITE(*,*)  JSXT_REAL 
	WRITE(*,*)  JSXT_IMAG 
      WRITE(*,*)  JSYT_REAL 
	WRITE(*,*)  JSYT_IMAG

	WRITE(*,*) 'CALCULATE THE GROUND CURRENT FOR ANOTHER PLANE: (Y/N)?'
	READ(*,*) ANSWER_AGAIN

		IF (ANSWER_AGAIN=='Y' .OR. ANSWER_AGAIN=='Y')THEN
			REPEAT = 1;
		ELSE
			REPEAT = 0;
		END IF

	END DO 

	END IF
	
	pause

	END PROGRAM MOM